%**************************************************************************
% DISE-2024 with CasADi
% Baseline calibration (No-intervention)
% A. Bongers and J. L. Torres (2025)
% Model: DISE-2024
% This version uses the IPOPT algorithm without terminal conditions
%**************************************************************************
clc; 
clear all;
close all;

% Model parameters
Params  =   Space1999;
Params.T    =   250;
Params.beta =   1/(1+Params.rho);
TT = 50;
% Solution
    sol =   ocp(Params);
    k   =   sol.k;
    k   =   k(1:Params.T-TT);
    s   =   sol.s;
    s   =   s(1:Params.T-TT);
    c   =   sol.c;
    c   =   c(1:Params.T-TT);
    is  =   sol.is;
    is  =   is(1:Params.T-TT);
    ik  =   sol.ik;
    ik  =   ik(1:Params.T-TT);
    W   =   sol.W;
    W   =   W(1:Params.T-TT);
    Z   =   sol.Z;
    Z   =   Z(1:Params.T-TT);
    F   =   sol.F;
    F   =   F(1:Params.T-TT);

% Auxiliary variables
    for i=1:Params.T-TT
        D(i) = W(i)+Z(i)+F(i);    
        y(i) = Params.a(i)*k(i)^Params.alpha1*s(i)^Params.alpha2*Params.N(i)^(1-Params.alpha1-Params.alpha2);
        L(i) = Params.mu*(1-Params.m(i))*is(i)/Params.eta;
        S(i) = Params.mu*s(i)/Params.q(i);
        X(i) = Params.theta*D(i)*S(i);
        KSS(i)=Params.theta*D(i);
    end

% Figures
Periods=2023:2023+Params.T-1-TT;
figure;
subplot(2,2,1); plot(Periods, k, 'b', 'LineWidth', 2);
xlabel('Time'); ylabel('Earth Capital'); title('Earth Capital Path'); grid on;
subplot(2,2,2); plot(Periods, s, 'b', 'LineWidth', 2);
xlabel('Time'); ylabel('Space Capital'); title('Space Capital Path'); grid on;
subplot(2,2,3); plot(Periods, c, 'r', 'LineWidth', 2);
xlabel('Time'); ylabel('Consumption'); title('Consumption Path'); grid on;
subplot(2,2,4); plot(Periods, y, 'r', 'LineWidth', 2);
xlabel('Time'); ylabel('Output'); title('Output Path'); grid on;

figure;
subplot(2,2,1); plot(Periods, L, 'b', 'LineWidth', 2);
xlabel('Time'); ylabel('Launches'); title('Launches'); grid on;
subplot(2,2,2); plot(Periods, S, 'b', 'LineWidth', 2);
xlabel('Time'); ylabel('Satellites'); title('Operational satellites'); grid on;
subplot(2,2,3); plot(Periods, X, 'r', 'LineWidth', 2);
xlabel('Time'); ylabel('Satellites'); title('Satellites destroyed'); grid on;
subplot(2,2,4); plot(Periods, KSS, 'r', 'LineWidth', 2);
xlabel('Time'); ylabel('Probability'); title('The Kessler syndrome'); grid on;

figure;
subplot(2,2,1); plot(Periods, D, 'b', 'LineWidth', 2);
xlabel('Time'); ylabel('Debris'); title('Orbital debris > 1 cm'); grid on;
subplot(2,2,2); plot(Periods, W, 'b', 'LineWidth', 2);
xlabel('Time'); ylabel('satellites'); title('Derelict satellites'); grid on;
subplot(2,2,3); plot(Periods, Z, 'r', 'LineWidth', 2);
xlabel('Time'); ylabel('Rocket bodies'); title('Rocket bodies'); grid on;
subplot(2,2,4); plot(Periods, F, 'r', 'LineWidth', 2);
xlabel('Time'); ylabel('Fragments'); title('Debris fragments'); grid on;
 